﻿/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Jurgen Van Gael
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
using System;
using dnAnalytics.Math;
using dnAnalytics.Properties;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Decomposition;

namespace dnAnalytics.Statistics.Distributions
{
    /// <summary>
    /// This class implements functionality for the inverse Wishart distribution. This distribution is
    /// parameterized by the degrees of freedom nu and the scale matrix S. The inverse Wishart distribution
    /// is the conjugate prior for the covariance matrix of a multivariate normal distribution.
    /// </summary>
    /// <remarks>The distribution will use the <see cref="System.Random"/> by default. 
    /// Users can set the random number generator by using the <see cref="RandomNumberGenerator"/> property.
    /// 
    /// The statistics classes will check all the incoming parameters whether they are in the allowed
    /// range. This might involve heavy computation. Optionally, by setting Control.CheckDistributionParameters
    /// to false, all parameter checks can be turned off.</remarks>
    public class InverseWishart
    {
        // The degrees of freedom for the inverse Wishart distribution.
        private readonly double mNu;

        // The scale matrix for the inverse Wishart distribution.
        private readonly Matrix mS;

        // Caches the Cholesky factorization of the scale matrix.
        private readonly Cholesky mChol;

        /// <summary>
        /// Construct a new inverse Wishart distribution.
        /// </summary>
        /// <param name="nu">The degrees of freedom for the inverse Wishart distribution.</param>
        /// <param name="S">The scale matrix for the inverse Wishart distribution.</param>
        public InverseWishart(double nu, Matrix S)
        {
            mNu = nu;
            mS = S;
            mChol = new Cholesky(mS);

            if (Control.CheckDistributionParameters)
            {
                CheckParameters(nu, mChol);
            }

            RandomNumberGenerator = new System.Random();
        }

        /// <summary>
        /// A string representation of the distribution.
        /// </summary>
        public override string ToString()
        {
            return "InverseWishart(Nu = " + mNu + ", Rows = " + mS.Rows + ", Columns = " + mS.Columns + ")";
        }

        /// <summary>
        /// Checks the parameters of the inverse Wishart distribution.
        /// </summary>
        /// <param name="nu">The degrees of freedom for the inverse Wishart distribution.</param>
        /// <param name="S">The scale matrix for the inverse Wishart distribution.</param>
        /// <exception cref="ArgumentOutOfRangeException">If the degrees of freedom parameter is negative.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If the scale parameter is not positive definite.</exception>
        public static void CheckParameters(double nu, Matrix S)
        {
            CheckParameters(nu, new Cholesky(S));
        }

        /// <summary>
        /// Checks the parameters of the inverse Wishart distribution.
        /// </summary>
        /// <param name="nu">The degrees of freedom for the inverse Wishart distribution.</param>
        /// <param name="chol">The cholesky factorization of the scale matrix for the inverse Wishart distribution.</param>
        /// <exception cref="ArgumentOutOfRangeException">If the degrees of freedom parameter is negative.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If the scale parameter is not positive definite.</exception>
        public static void CheckParameters(double nu, Cholesky chol)
        {
            if (nu < 0.0)
            {
                throw new ArgumentOutOfRangeException("nu", Resources.NotPositive);
            }

            // Try a matrix factorization of S.
            if (!chol.IsPositiveDefinite())
            {
                throw new NotPositiveDefiniteException();
            }
        }

        /// <summary>
        /// The degrees of freedom for the inverse Wishart distribution.
        /// </summary>
        public double Nu
        {
            get { return mNu; }
        }

        /// <summary>
        /// The scale matrix for the inverse Wishart distribution.
        /// </summary>
        public Matrix S
        {
            get { return mS; }
        }

        /// <summary>
        /// Gets the mean.
        /// </summary>
        /// <value>The mean.</value>
        public Matrix Mean
        {
            get { return mS/(mNu - mS.Rows - 1.0); }
        }

        /// <summary>
        /// Samples an inverse Wishart distributed random variable by sampling
        /// a Wishart random variable and inverting the matrix.
        /// </summary>
        public Matrix Sample()
        {
            return Sample(RandomNumberGenerator, mNu, mS);
        }

        /// <summary>
        /// Samples an inverse Wishart distributed random variable by sampling
        /// a Wishart random variable and inverting the matrix.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="nu">The degrees of freedom.</param>
        /// <param name="S">The scale matrix.</param>
        public static Matrix Sample(System.Random rnd, double nu, Matrix S)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(nu, S);
            }

            Matrix R = Wishart.Sample(rnd, nu, S.Inverse());
            return R.Inverse();
        }

        /// <summary>
        /// Gets or sets the random number generator.
        /// </summary>
        /// <value>The random number generator used to generate a random sample.</value>
        public System.Random RandomNumberGenerator { get; set; }


        /// <summary>
        /// Evaluates the probability density function for the inverse Wishart distribution.
        /// </summary>
        /// <param name="X">The matrix at which to evaluate the density at.</param>
        /// <exception cref="ArgumentOutOfRangeException">If the argument does not have the same dimensions as the scale matrix.</exception>
        public double Density(Matrix X)
        {
            int p = mS.Rows;

            if (X.Rows != p || X.Columns != p)
            {
                throw new ArgumentOutOfRangeException("X", "X must be of same size as the scale matrix.");
            }

            Cholesky chol = new Cholesky(X);
            double dX = chol.Determinant();
            Matrix SXi = chol.Solve(S);

            // Compute the multivariate Gamma function.
            double gp = System.Math.Pow(Math.Constants.Pi, p*(p - 1.0)/4.0);
            for (int j = 1; j <= p; j++)
            {
                gp *= SpecialFunctions.Gamma((mNu + 1.0 - j)/2.0);
            }

            return System.Math.Pow(dX, -(mNu + p + 1.0)/2.0)
                   *System.Math.Exp(-0.5*SXi.Trace())
                   *System.Math.Pow(mChol.Determinant(), mNu/2.0)
                   /System.Math.Pow(2.0, mNu*p/2.0)
                   /gp;
        }
    }
}